#include <stdio.h>
#include <assert.h>
#include <math.h>
#define N 1000

double f1(double x,double a,double b)
{
	assert(x >= a && x <= b);
	return 1/x;
}

double f2(double x,double a,double b)
{
	assert(x >= a && x <= b);
	return 4/(1+x*x);
}

double Simpson(double a,double b,int n,double (*f)(double ,double ,double ))
{
	double h = 0.001;
	double x[N+1];
	double s = 0;
	int i;
	for(i = 0;i < n+1;i++)
		x[i] = a + i*h;
	s += f(a,a,b);
	s += f(b,a,b);
	for(i = 0;i < n;i++)
		s += 4*f((x[i]+x[i+1])/2,a,b);
	for(i = 1;i < n;i++)
		s += 2*f(x[i],a,b);
	s = s*h/6.0;
	return s;
}

int main()
{
	double s;
	double t;
	s = Simpson(1,2,N,f1);
	t = log(2.0);
	printf("Simpson ln2：%.16f\n",s);
	printf("Math ln2：%.16f\n",t);
	s = Simpson(0,1,N,f2);
	t = 3.141593;
	printf("Simpson PI：%.16f\n",s);
	printf("Math PI：%.16f\n",t);
	return 0;
}
